Critical line of honeycomb-lattice anisotropic Ising antiferromagnets in a field 
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We use numerical transfer-matrix methods, together with finite-size scaling and conformal in- 
variance concepts, to discuss critical properties of two-dimensional honeycomb-lattice Ising spin- 1/2 
magnets, with couplings which are antiferromagnetic along at least one lattice axis, in a uniform 
external field. We focus mainly on the shape of the phase diagram in field-temperature parameter 
space; in order to do so, both the order and universality class of the underlying phase transition are 
examined. Our results indicate that, in one particular case studied, the critical line has a horizontal 
section (i.e. at constant field) of finite length, starting at the zero-temperature end of the phase 
boundary. Other than that, we find no evidence of unusual behavior, at variance with the reentrant 
features predicted in earlier studies. 

PACS numbers: 64.60.De, 75.10.Hk, 75.30.Kz 



I. INTRODUCTION 

The work reported in this paper is a study of Ising 
spin-1/2 systems on a honeycomb (HC) lattice, with 
ferro- (F) and antiferromagnetic (AF) interactions, in 
the presence of a uniform magnetic field. With k = 1, 2, 
and 3 being the three lattice directions, the Hamiltonian 
is given by: 
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where the subscript (i, j) k denotes nearest-neighbor spins 
along lattice direction k, and <Jij = ±1. Here, all fields 
H, coupling strengths J^, and temperatures T will be 
given in units of J± . We only consider cases for which at 
least one of the Jk is negative, thus the AF character is 
always present. 

It is well known that the two-dimensional Ising model 
is not amenable to exact solution when an external field 
is introduced. If the field is uniform and the system is 
AF, an ordered phase is found for suitably low values 
of H and T. When the zero-field ground state exhibits 
macroscopic entropy, the phase diagram on the H — T 
plane shows reentrant behavior. This happens for the 
triangular-lattice isotropic AF [l|, for which it has been 
found that the second-order transition along the critical 
line belongs to the ferromagnetic three-state Potts uni- 
versality class 

If the H = ground state has vanishing residual en- 
tropy per spin, a simpler picture is expected to hold, in 
which the Ising character of the zero-field transition is 
preserved everywhere along the critical line on the H — T 
plane. It has been shown to be so, for isotropic AFs on 
both the square and HC (i-tToj lattices. 

Anisotropic square-lattice Ising systems with mixed in- 
teractions (F along one lattice direction, AF along the 
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other) in a field have been discussed extensively; see 
Ref. [ll| for a summary of early work, and Ref. [l2| for 
more recent references. In this case, similarly to the pure 
AF the ground-state residual entropy vanishes; remark- 
ably, reentrant behavior was found at low temperatures 
in some numerical or analytic treatments. While a non- 
trivial ground-state structure is not known to be a pre- 
requisite for the existence of reentrances, the latter type 
of phenomenon usually results from subtle combinations 
of competing effects, whose nature it would be of impor- 
tance to establish. The results of Ref. indicate that the 
critical line starts horizontally at the zero-temperature 
end of the phase boundary, thus excluding reentrant be- 
havior for this system. 

For HC lattices with anisotropy, the existence of reen- 
trances has been predicted [l3| for a variety of combina- 
tions of F and AF interactions, depending also on their 
relative strength. In Ref. Il3l . an approach was used which 
considers the zeros of the partition function on an ele- 
mentary lattice cycle, and their connection to the free 
energy singularity at the transition §. When applied 
to the mixed square-lattice model described in the pre- 
ceding paragraph, it also predicted a reentrant critical 
curve |13 |. 

Our purpose here is to establish, by numerical meth- 
ods, the shape of the phase diagrams of systems described 
by Eq. fT]), especially as regards the existence (or not) of 
reentrant sections; in order to do so, we necessarily probe 
the universality class of the phase transitions along the 
respective critical lines. 

We use transfer-matrix (TM) methods, in connection 
with finite-size scaling and conformal invariance ideas, in 
order to produce numerically accurate phase diagrams. 
The underlying hypotheses in our work are: (i) that the 
phase transition is second-order all along the critical line, 
and (ii) that it belongs to the Ising universality class. 
Both assumptions are critically reviewed toward the end 
of the paper, in light of the numerical results obtained 
while assuming their validity. 

In Sec. [H] we recall the calculational methods used for 
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the approximate location of the critical line. In Sec. IIIII 
the respective results are exhibited. In Sec. IIV1 we ana- 
lyze the data generated in Sec. IIIII both in comparison 
with the existing literature, and in regard to their inter- 
nal consistency; finally, concluding remarks are made. 



II. CALCULATIONAL METHOD 
A. Introduction 

Referring to Eq. |T]) , here we consider combinations of 
interaction signs in which either one, two or all three of 
the Jk are AF. The respective strengths will be specified 
below in each case, in order to reproduce points in {Jk} 
parameter space for which Ref. [l3| predicts sizable reen- 
trant sections of the respective critical line. As explained 
in the following, it is possible to do so while always keep- 
ing the couplings along two directions with the same sign 
and strength, while bonds along the third direction (to 
be denoted as inhomogeneous) differ from the other two 
in strength and/or sign. 

We set up the TM on strips of width N sites, with 
periodic boundary conditions across. Here we use two 
distinct choices of orientation, respective to the lattice 
axes: in (a) the TM proceeds perpendicularly to one lat- 
tice direction |14| . while in (b) it goes parallel to one 
lattice direction For (a), periodicity across the strip 
imposes that N must be even; no parity restrictions arise 
for (b). We kept to the range 4 < N < 20, which (to- 
gether with suitable extrapolation techniques) generally 
proved enough to yield accurate estimates of the critical 
lines. For the present case of anisotropic systems, and 
bearing in mind the definition of inhomogeneous bonds 
introduced in the preceding paragraph, it is possible to 
consider the following variants: 

(al) choice (a), with the inhomogeneous bond perpendic- 
ular to the TM's direction of advance; 
(a2) (a), with the inhomogeneous bond along either of 
the remaining two directions; 

(bl) choice (b), with the inhomogeneous bond parallel to 
the TM's direction of advance; 

(b2) (b), with the inhomogeneous bond along either of 
the remaining two directions. 

As these are weakly anisotropic systems (TBI . Utij . one 
would expect estimates of, e.g., critical exponents and 
locations of critical points, to converge to the same 
orientation-independent limit for N 3> 1, while finite- 
size corrections should differ in each case. However, the 
essentially one-dimensional character of the strips used in 
our calculations, coupled with the distinct nature of spin 
couplings along the lattice axes, may introduce subtle bi- 
ases when one iterates the TM along a direction which is 
fixed with respect to those axes. Thus, it is important to 
consider alternatives which at least partially compensate 
for such possible distortions. We used two distinct pro- 
cedures, presented respectively in Sections III Bl and III CI 



of which the latter is expected to be less prone to the 
aforementioned biases than the former. 



B. Keeping rj = 1/4 

Following earlier work on similar problems [H], 0, 
UH, our finite- N estimates for the critical line are found 
by requiring that the amplitude-exponent relation of con- 
formal invariance on strips [I?} be satisfied, with the Ising 
decay-of-correlations exponent r\ = 1/4: 

4Nk n (T,H) = (tt , (2) 

where k^(T,H) = In | A1/A2 | is the inverse correlation 
length on a strip of width N sites, with Ai, A2 being the 
two largest eigenvalues (in absolute value) of the TM, 
and £ is a geometric factor which compensates for the 
fact that on the HC lattice the strip width, in lattice 
parameter units, is not equal to N. One has £ = 2/\/3 
for orientation (a) , and £ = 1/ \/3 for (b) . 

C. Phenomenological Renormalization 

The assumption made in Eq ©, that the phase tran- 
sition belongs to the Ising universality class, can be re- 
laxed. Instead, one can demand only that it remain of 
second order. From finite-size scaling, one gets the basic 
equation of the phenomenological renormalization group 
(PRC) [H| for the critical line: 

Nk n (T,H) = N'k n ,(T,H), (3) 

where the strip widths N and N' are to be taken as close 
as possible for improved convergence of results against 
increasing TV. 

Note that for PRG one is always comparing correla- 
tion lengths evaluated along the same lattice direction 
in Eq. ©, so the likely biases mentioned above tend to 
cancel out, if present |l2j. 

PRG results can also be used as a test of the inter- 
nal consistency of the Ising- universality class assumption; 
furthermore, should additional, non-Ising, transitions be 
present elsewhere in parameter space, they should be de- 
tected by PRG. 

III. RESULTS 

A. Ji = J 2 < 0; J 3 > 

Analysis of zero-field ground-state configurations 
shows that, in this case, N must be restricted to mul- 
tiples of 4 for choice (a2) above, and to even values for 
choice (b2). 

We take J3 = 1. Consequently: (i) ground-state con- 
siderations show that H = 2 is the zero-T critical field, 
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Figure 1. (Color online) For Ji = J2 = — J3 = —1, low- 
temperature approximate critical boundaries given by solu- 
tions of Eq. ([2]) [empty symbols: TV = 4; full symbols, TV = 14 
(al, bl), 12 (a2,b2)], or Eq. © [PRG_al: strips of widths 
TV, TV — 2, TV = 12, geometry al]. Large- TV curves for a2 and 
bl practically coincide on the scale shown. See text. 



and (ii) the zero-field critical temperature coincides with 
that of the pure AF (or F) system: T C (H = 0) = 
2/ln(2 + \/3) = 1.5186514.... 

In Ref. [l3| it is predicted that, for J2/J1 = | J3/J1I > 
0.6725 the critical curve starts from H = going to- 
wards higher temperatures (thus forming a "bulge") and 
then turns towards lower T, monotonically approaching 
its limiting intercept at H = H C (T = 0). Further- 
more, the critical curve is predicted to reach the point 
T = 0, H = H C (T = 0) horizontally. 

By numerically solving Eqs. (J2JI and (J3J), we found that 
for H <C 1, the approximate critical curves leave the T 
axis vertically, and are very close to each other, for all 
possible choices (al)— (b2) described in Section [IT] We 
find no evidence of a "bulge": T C (H) decreases monoton- 
ically with increasing H. We shall return to this point 
below. As H increases, differences between curves gen- 
erated by the various procedures become slightly more 
pronounced. These are illustrated in Fig. [T] One sees 
that, with increasing TV, all four families of solutions of 
Eq. Q converge towards approximately the same inter- 
mediate location, which also coincides with the single 
solution of Eq. ((3]) shown. We have found that the solu- 
tions of Eq. © converge much faster with increasing TV 
than those of Eq. @, so the single PRG curve depicted 
accurately represents the TV — > 00 limit of its family, to 
the scale of the figure. 

At low temperatures T < 0.2 — 0.3, numerical difficul- 
ties arise, because of the very large ratio between posi- 



Figure 2. (Color online) Approximate phase diagrams gener- 
ated by the solutions of Eq. ©, in geometry (al) for N=16, 
for Ji — J2 ~ —J3 = —1 (squares), and for 2Ji = — J2 = 
— J3 = 2 (hexagons). 



tive and negative exponentials which are the TM states' 
Boltzmann weights in that region. Even accounting for 
this, one can see that all approximate critical curves un- 
equivocally home in toward (T, H) = (0, 2) at finite an- 
gles. This is illustrated by the dashed straight-line seg- 
ments in Figure [T] which are guides to the eye only but 
most likely are good stand-ins for the actual shape of 
the corresponding curves. In conclusion, we do not find 
any evidence that the critical curves approach the T = 
horizontally. 

Going back to the predicted "bulge "-like behavior pre- 
dicted at higher temperatures, Fig. [2] shows the full 
phase diagram generated by the solutions of Eq. ©, 
in geometry (al) for TV = 16, both for the case pre- 
viously discussed of J\ = J 2 = —J3 = — 1, and also 
for 2Ji = — J2 = —J3 = 2. This second combination 
has been considered because it allows direct comparison 
with one of the diagrams shown in Figure 12 of Ref. Il3l . 
corresponding to the same coupling values. The phase 
diagram obtained there exhibits a horizontal portion at 
H = 4, from T = to T w 1, as well as a "bulge" with 
maximum extent at T « 2.7, H ps 2.3. Here, both fea- 
tures are absent from the numerically evaluated critical 
lines. 



B. Ji = J 2 > 0; J 3 < 

As in the case of Section MI Al here TV must be a mul- 
tiple of 4 for choice (a2) above, and even for choice (b2). 
We assume J3 = —1, so (i) the zero-temperature crit- 
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Figure 3. (Color online) For Ji = J2 ~ —J3 = 1, low- 
temperature approximate critical boundaries given by solu- 
tions of Eq. ([2]) [empty symbols: TV — 4; full symbols, TV = 14 
(al), 16 (a2,bl), 12 (b2)], or Eq. © [PRGal: strips of 
widths N, N - 2, N = 12, geometry al]. See text. 



Figure 4. (Color online) For Ji — J2 = — J3 = 1, excess 
peak heights (above H = 1) of numerically-obtained phase 
diagrams from solutions of Eq. ([2]), for geometries (al), (a2), 
and (b2), against inverse strip width 1/JV. See Figure [3] for 
reference. 



ical field is H C (T = 0) = 1, and (ii) T C {H = 0) = 
2/ln(2 + V3) as in Section HTlAl 

In Ref. |l3| it is predicted that, for J2/J1 = \Ja/Ji\ > 
1/3, the critical line should leave the T — axis with 
positive slope. For J3 = —1, from Figure 11 of Ref. [TH 
one sees that the peak of the corresponding reentrance is 
expected to occur at T m 0.45, H w 1.05. No "bulge", i.e. 
a section of the critical curve extending to T > T C {H = 0) 
at low H, is predicted. 

Our results for low T, encompassing the region of the 
predicted reentrance, are shown in Fig. [3] For geome- 
tries (al), (a2), and (b2), our numerical results indeed 
show reentrant-like behavior in the predicted range of T. 
However, in all three cases the evolution of the excess 
peak heights (i.e. above the H = 1 level) against in- 
creasing N is towards reduction. This is illustrated in 
Figure [3J It is clearly seen that the trend followed in all 
cases certainly excludes a positive limiting height; on the 
contrary, slightly negative values (sa 0.005 in modulus) 
appear more likely. 

The solutions of Eq. @ in geometry (bl) do not 
show reentrances; their low-temperature sections ap- 
proach straight lines homing in towards (T,H) = (0, 1). 
Upon increasing N, the negative slope of such straight- 
line sections of the T c x H curves becomes smaller in ab- 
solute value. A similar trend is followed by the solutions 
of Eq. ((3]) in the same geometry. 

Both sets of slope values are exhibited in Figure [5] For 
the solutions of Eq. ([2]) a simple linear extrapolation of 
large-iV data gives a negative limiting slope w —0.030. 




1/N 



Figure 5. (Color online) For Ji = Jg = — J3 = 1, slopes of the 
low-temperature straight-line sections of approximate critical 
boundaries given in geometry (bl), by the solutions of Eq. © 
or Eq. (|3j) [ PRC ], against inverse strip width. Data points are 
for N = 6 - 20. Lines for TV -1 < 0.05 are fits of 15 < N < 20 
data to forms shown [ "par." stands for "parabolic" ]. 
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Attempts to include a quadratic term (not shown) result 
in a limiting value around —0.015. Fitting to a single- 
power form ~ TV~ X , thus forcing the limiting slope to 
vanish as TV — > oo, requires x w 0.69. Though not alto- 
gether implausible, an exponent x < 1 means a growing 
amount of curvature with increasing TV. Taken together, 
the preceding considerations indicate that a positive ini- 
tial slope of the critical curve appears unlikely. The solu- 
tions of Eq.® in geometry (bl) behave smoothly, and a 
parabolic fit [i.e., with both linear and quadratic terms 
in 1/TV] gives a limiting slope equal to (6 ± 2) x 10~ 4 , 
which essentially equates to zero in the present context. 

On the other hand, PRG estimates in geometry (al) 
[shown in Figure [3]], and also in (a2) and (b2) [not 
shown ] consistently give the approximate critical bound- 
ary lying slighty below the H = 1 line (within less than 
0.1% of it) for all T < 0.4. As remarked in Section ImAl 
here too the solutions of Eq. ^ [ other than those for ge- 
ometry (bl) ] exhibit very little TV- dependence: for (al), 
differences between results for TV — 6 and TV — 12 are at 
most of order 1 — 2 parts in 10 3 , with the largest values oc- 
curring midway between the phase diagram's endpoints. 
Thus, one cannot discard the possibility that the criti- 
cal curve starts horizontally from (T,H) = (0,1), and 
remains flat for a finite extent, up to T ss 0.4. 
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Figure 6. (Color online) For Ji — J2 = —I, J'i = —0.4, 
examples of approximate critical boundaries from solutions of 
Eq. © (all for TV = 8), and Eq. © (with TV = 10, TV' = 8). 



C. Ji = J 2 < 0; J 3 < 

In this case, the zero-field ground state coincides with 
that of the pure AF, thus there are no additional re- 
strictions on strip widths other than those mentioned in 
Section UTAl 

We assume J3 = —0.4, so (i) the zero-temperature crit- 
ical field is H C {T = 0) = 2.4, and (ii) T c ° = T C {H = 0) 
does not take on the exact value 2/ ln(2 + as in Sec- 
tions hie] and mm 

For these coupling values, Ref. [l3| predicts that the 
critical curve should leave the T = axis with a posi- 
tive slope, S = § In 2, see their Eq. (38) and following 
considerations. 

The eight sets of finite-N estimates, obtained by solv- 
ing Eqs. @ and ([3]) in all four geometries (al)— (b2), 
give very similar results, as depicted in Figure [6l In 
contrast with Sections MI Al and MI Bl the largest vari- 
ations among differing calculational schemes are found 
for H < 1. At H = we quote T c ° = 1.1170(5), where 
the error bar reflects the scatter among extrapolations of 
finite- TV sequences for each of the eight sets available. At 
T = the critical curve starts with a negative slope, in 
disagreement with the prediction of Ref. Il3l . 



IV. DISCUSSION AND CONCLUSIONS 

For each of the three specific combinations of coupling 
signs and strengths discussed in Section Mil we generally 



found sets of solutions to Eqs. @ and © in good agree- 
ment with each other, and clearly showing convergence 
towards the same asymptotic location of the respective 
phase boundaries for increasing TV. We only failed to 
find such solutions for low temperatures, T < 0.3 in the 
worst cases. However, even in that limit we always man- 
aged to get enough data to show that the approximate 
critical curves are homing in towards the known zero- 
temperature critical field. Taken together, these consid- 
erations indicate that the phase transition between or- 
dered and paramagnetic states remains second-order ev- 
erywhere along the phase boundary [thus justifying the 
use of Eq. (J3|) ] , and that the transition is in the Ising 
universality class [ which supports Eq. @ ] . 

With the possible exception of the low-temperature re- 
gion in Section MI Bl (to be discussed separately below) , 
we found no evidence of reentrances, "bulges", or hori- 
zontal sections in the numerically-calculated phase dia- 
grams here presented. Near T = our curves leave the 
H axis with (negative) slope S, as defined in: 

H c (T) = H c (0)+ST , (4) 

and for H — > their shape is well fitted by a parabolic 
form: 

T c {H)=T c (0)-aH 2 . (5) 

Our estimates of S and a are given in Table fl] For each 
set of couplings, the error bars reflect the scatter among 
fits of large-TV approximate curves, each corresponding 
to one the eight possible combinations of calculational 
procedure and geometry. The comparison of calculated 
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Table I. Adjusted values of initial slope S at T — [followed 
by predictions from Ref. HU], and quadratic term a in 
parabolic fit at H — 0, for phase diagrams in Sections IIII Af - 
IIII CI [see, respectively, Eqs. © and ©]. 



Section 


S 


S (Ref. 13) 


a 


|IIIA| 


-0.321(7) 





0.205(3) 


|IIIB| 


£ [-0.03, +6 x 10" 4 ] 


(1/4) In 2 


0.504(2) 


Unci 


-0.995(5) 


(3/2) In 2 


0.158(1) 



initial slopes with the predictions of Ref. 8, summarized 
in Table HJ bears strong resemblance to the discussion of 
Ref. UH for similar sytems on a square lattice. Also, as re- 
marked there, it is worth recalling the comparable prob- 
lem of isotropic antiferromagnets (on both square and HC 
lattices): although the critical lines H C {T) given in Ref. 8: 
do not exhibit reentrances, they are always above those 
found in Refs. and (except at the T — and H = 
ends, where the lines coincide in both cases). Thus the 
results presented here are consistent with previous ones, 
in indicating that the methods employed in Refs. HI and 
Il3l tend to overestimate the extent of the ordered region 
in parameter space. 

Going back to the low-temperature region of the phase 
diagram in Section llll A[ the behavior seen there strongly 
suggests that the critical line approaches the T — axis 
at a very low angle, possibly even horizontally or (though 
less likely) at a very slight reentrance. The range of slope 
estimates given in Table U reflects the arguments given 



in connection with Figures [3J HI an d [5] Similar behav- 
ior occurs in the square-lattice metamagnet studied in 
Ref. [Ill, where a sizable extent of the low-temperature 
phase boundary is found to be horizontal to within 0.1% 
(the same fractional deviation exhibited by the PRG 
curves here, see Section fill Bl) . 

In the three-dimensional version of the problem of 
Ref. [TU, it is known that a tricritical point occurs at low 
temperatures [l^. As shown in Figure [3J we managed to 
find solutions of Eq. ([3]) [whose validity is based on the 
assumption of a second-order transition] for rather low 
T . Furthermore, in that region the approximate critical 
curves are well-behaved as far down in temperature as 
we can ascertain, pointing directly towards the T — 
critical field. Thus no evidence has been found that the 
nature of the transition changes, or is about to change, 
close to T = 0. 

In conclusion, we have found strong indications that 
reentrant behavior does not occur in the systems studied 
here, for which the ground state has a relatively simple 
structure. On the other hand, even for these systems 
it is possible to find slightly unusual features such as 
horizontal sections of the phase diagram, which it do not 
obviously suggest themselves. 
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